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SUMMARY 

The three-dimensional compressible Navler Stokes equations are formulated 
In a '•otatlng coordinate system, so as to Include centrifugal and Coriolis 
forces. The equations are parabolized by using a previously calculated Invls- 
c 1 d static pressure field. The thin layer Navler-Stokes approximation, which 
neglects streamwlse diffusion, Is used. A body-fitted coordinate system Is 
used. The streamwlse momentum equation Is uncoupled from the cross-stream 
momentum equation by using contravarlant momentum components, and then using 
the contravar lant velocity components as primary unknowns. To reduce problems 
with small separated regions, the Reyhner and rlugge-Lotz approximation Is 
used. The energy equation Is Included to allow for calculation of heat trans- 
fer. The flow may be laminar, or a simple eddy- vl scos 1 ty turbulence may be 
used. A number of curved ducts and an axial stator have been analyzed, 
Including cases for which experimental data are available. 
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INTRODUCTION 

Well-guided Internal flow without separated regions behaves In a para- 
bolic manner (ref. 1) so that the flow field can be predicted using a para- 
bolic, streamwlse marching method. Several computer codes have been written 
for calculating three-dimensional viscous flow through ducts of fairly complex 
geometry (rr c '. 2 and 3). The same principles can be applied to calculate 
flow through the blades of a turbomachine. However, there are several distin- 
guishing aspects of flow through a turbomachine blade row: (1) The blade row 
may be rotating. (2) The hub and shroud are usually surfaces of revolution. 

(3) The passage has four walls meeting at an angle. (4) The flow may be peri- 
odic upstream and downstream of the blade. (5) The shroud Is usually station- 
ary for a rotating blade row. The method and computer code described here Is 
especially formulated to satisfy and take advantage of these particular 
aspects of the flow. 
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The parabolic marching calculation Is analagcus to calculation of a 
three-dimensional boundary layer. An Invlscld pressure field Is used as an 
Initial pressure field; then a single-pass calculation through the length of 
the blade passage Is made. To simplify the application of boundary conditions 
on complex blade surfaces, a body-fitted coordinate system Is used. This 
body-fitted coordinate system has one coordinate which generates circles about 
the axis of rotation. One of the coordinate surfaces coincides with the hub of 
the blade row, and another coincides with the shroud. This type of coordinate 
system eases Implementation of periodic boundary conditions upstream and down- 
stream of the blade row; on the other hand, this type of coordinate system 
cannot be orthogonal. 

In the spirit of Patankar and Spalding (ref. 1) the streamwlse momentum 
calculation Is uncoupled from the cross-stream momentum calculation. This 
uncoupling Is accomplished by using contravarlant components of the momentum 
equation, and using contravarlant velocity components as the primary unknowns. 
The downstream contravarlant velocity components are always calculated Impli- 
citly for stability. To calculate the streamwlse contravarlant velocity com- 
ponents at the next downstream station, It Is assumed that the static pressure 
and cross- stream contravarlant velocity components are known at the previous 
upstream station. The streamwlse pressure gradient Is varied by a constant 
amount over the passage cross-section ro that global mass flow through the 
entire passage Is conserved. The cross-stream contravarlant velocity compo- 
nents are calculated at the next downstream station based on the Invlscld 
pressure field, and using the last calculated streamwlse contravarlant veloc- 
ity components. Continuity Is then checked for each mesh region, and the 
cross- stream pressure gradient Is adjusted so that continuity Is satisfied for 
every mesh region. Iteration Is required at each station so that both cross- 
stream momentum and local continuity are satisfied. For turbulent flow, the 
Baldwin- Lomax eddy - vl scos Ity turbulence model Is used (ref. 4). 


PARA80I IC nARCHING THROUGH A TURBOMACHINE PASSAGE 

Contravarlant -'orm of the Navler- Stokes Equations for a Rotating 
Non- orthogonal Coordinate System 

The continuity equation In conservation form In cylindrical coordinates 
Is: 

ft < r "> ' fr < r 'V * fe '“V "t; < r »«z> ■ 0 <’> 

In cylindrical (r,e,z) coordinates the Navler- Stokes equations may be written 
In conservation form as follows: 

ft C'V * fr * fi <-> V rV ‘ fe (^W 

pV e " 5r ^ rT rr “ rp ^ * Je ^ T re^ * az ^ rT rz^ T ee * p 
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' pV r V 8 " If (r, re> * fe (, ee p > * If <r 'ez> * Ve 
If <» v z> ♦ If ^VV f ie < p Vz> - If ( p ev|) 
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The change to a rotating coordinate system Is accomplished by the change of 
variables : 


W r = V r 
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The partial derivatives with respect to t then become: 

3 


a_ a 

at = at' 


(•> 


30 


The resulting Navler- Stokes equations for rotating cylindrical coordinates are 
now {the ' has been dropped): 

It (rpw r» * If < rpw f> * fe ‘-’Ve 1 * If 0“rV 

- '•“I • r if - ’ee * If (r Vr> ‘ If <Ve> * If (r Vz> (2 > 
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( 3 ) 
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Finally, the energy equation In conservation form Is: 

+ V*^(e + p) =* v»k 7 T 7 * ( t 9 ) 

This Is expanded In cylindrical coordinates, and transformed tc rotating 
coordinates using 

a_ a a_ 

at “ at 1 '“ae 

Also, rothalpy I Is Introduced, using the relation 

e + p = pf ♦ wrpV 

The energy equation In conservation form, for rotating cylindrical coordinates 
and In terms of rothalpy Is: 

1^; (r(pl - p)j + r7* ( plfi) = rv*kVT <• r7*(t fl) 

or 
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Equations (1) to (5) can be written In compact form as 

a^q *■ 3pE + 3 qF + a z G = K + a r R + 3 qS + a z T (6) 
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Equation (6) can be transformed to an arbitrary nonorthogona 1 rotating 
coordinate system (E, n, f) using the chain rule. When this Is done, consider- 
able simplification occurs by combining terms properly (refs. 5 and 6), and by 
using the contravarlant velocity components as variables. Note that terms 
such as 
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are equal to zero (ref. 7), so that factors like 5 r / J 2 can be moved Inside 
the derivative. The form of the equations Is now: 
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Equation 7 still has cylindrical momentum components. It Is desired to 
utilize momentum components aligned with the mesh, so that streamwlse momentum 
Is easily uncoupled from the cross-stream momentum components. This Is accom- 
plished by taking the contravarlant momentum components. For example, the 
contravarlant ^-momentum component Is E r ( r-momentum) plus ri 0 /r( G-momentum) 
plus C z ( z-momentum) . This results In: 

ft* 8 ***** 

3 t q *■ 3j:E * a n F + a^G = a^R f a^S ^ 3 fT + K (8) 
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Derivation of the Navler-Stokes Equations for Parabolic Marching 
Through a Turbomachinery Blade Row 

Equation (8) retains alT terms from the original Navi er-Stokes , continu- 
ity, and energy equations for a general non orthogonal rotating coordinate 
system. Some terms will now be eliminated, or neglected. Only steady rela- 
tive flow Is considered so the time derivatives will be eliminated. Also, the 
thin layer assumption Is made so that all streamwlse derivatives of viscous 
terms will be neglected. Further, the mesh that Is used Is not completely 
general, so that certain coordinate derivatives will be eliminated. 

The coordinate system used Is shown In figure 1. The ^-coordinate Is In 
the streamwlse direction, n Is the blade to blade direction, and f Is from 
hub- to-shroud . Since the hub and shroud are both usually surfaces of revolu- 
tion, and to make It easier to apply periodic boundary conditions, the 
n-coordlnate lines will be circular arcs coincident with the e-coordlnate 
lines. Hence, £ and f are functions of r and z only, and sre Indepen- 
dent of e: 


1 1 


E - l(r.z), 


5e - 0 


n « n(r,e,z) 

C » C'r.z) 


The required coordinate derivatives, £ r , nr- ‘'tc. are calculated from 
the coordinate derivatives of tne Inverse trar.sf orm.' . *Om, r^, e^, etc., by the 
equations : 
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The Jacobian, J, Is given by: 
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Since 5 Is In the streamwlse direction, all E-derl vatl ves of the viscous 
terms are neglected. Further, by using the continuity equation, tor any 
function (< >, we have: 



Excluding continuity, the equations become: 
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Consistent with the thin boundary layer assumption, order of magnitude 
arguments are used to simplify calculation of the stress tensor In terms of 
the contravarlant velocity components. Terms that are Important Involve the 
second derivative of a velocity component parallel to a wall with respect to 
the distance from the sar »2 wall, whereas the derivative of a velocity compo- 
nent normal to a wall with respect to the same wall Is negligible. Also, 
coordinate derivatives (s r , n r * etc.) are considered to vary slowly compared 
to the velocity components themselves. Even with these simplifications, the 
calculation of the stress tensor Is rather complicated and messy for the non- 
orthogonal grid used. 


Numerical Solution Procedure 


A FORTRAN program called PARAMAR has been written to solve equation (12) 
by parabolic marching. Equation (12) Is the thin layer Navler-Stokes equation 
for a rotating, body fitted coordinate system. The ^-coordinate Is In the 
streamwlse direction, n Is blade-to-blade,and f Is hub-to-siiroud , as 
Indicated In figure 1. The equations are parabolized by using a previously 
calculated 1r;1sc1d static pressure field, and neglecting streamwlse diffusion 
terms. The pressure field may be obtained by the MERIDL and TSONIC codes 
(refs. J and 9) or from the Denton code (ref. 10). 

The form of the three momentum equations and the energy equation Is 


W* W* 1 3u> 3<p 

p r + p r an * p r ac 



S 


where may be any contravarlant velocity component or rothalpy. S Is 
the remaining terms of the equation. The coefficients a and b, and the 
source term S are different for each equation. By using a finite-difference 
’oproxlmatl on for the S-derlvatl ve, the equation may be approximated as: 


C, 
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It Is assumed that the mesh lines are coordinate lines, and that ?^ + -| - ^ * 1. 
Similarly, for the other coordinates, It will be assumed that An -= AC = 1- 
This equation Is an Implicit equation for In finite difference form, 

the equation may be solved for cpi+i directly by using the block trldlagonal 
algorithm (ref. 11, p. 196) 

Because of the Importance of accurately calculating continuity for Inter- 
nal flow, a staggered mesh Is used, much as suggested by Patankar and Spalding 
(ref. 1). The mesh arrangement and numbering sheme Is shown In figure 2 
Each of the momentum equations are uncoupled by using coefficients calculated 
at the upstream station. Since an 1 nvl sc 1 d pressure field Is used, correc- 
tions must be made to account for viscous effects In the boundary layer. 

The streamwlse (5) momentum equation Is used to calculate the stream- 
wise contravarlant velocity components at the next streamwlse station, using 
zero velocity boundary conditions. After calculating the downstream w£, the 
mass flow at the downstream station will be low due to an Increasing displace- 
ment thickness of the boundary layer. !ience, the pressure must be reduced at 
the downstream station, so that global continuity Is satisfied. This Is done 
by reducing the stream-wise pressure gradient uniformly over the entire pas- 
sage cross section. 

The cross-stream momentum equations are used to calculate the cross- 
stream contravarlant velocity components at the next station. Zero velocity 
boundary conditions are used. After calculating the cross-stream velocities, 
continuity Is checked for each mesh region. In general there will not be con- 
servation of mass for each mesh region. An equation can be derived to calcu- 
late a correction to the pressure gradient needed In the two cross-stream 
momentum equations to give cross-stream velocity components that will satisfy 
continuity for each mesh region. The resulting equation resembles a finite- 
difference Poisson equation for pressure over the downstream cross section. 

This correction to the pressure gradient Is needed to develop the correct 
cross flow within the boundary layer. Iteration Is required at each station to 
satisfy both the continuity equation and the momentum equations. 

If heat transfer Is negllble, constant rothalpy may be assumed. The 
energy equation Is used when heat transfer by convection and conduction Is of 
Interest. Then the rothalpy at the downstream station Is calculated Implic- 
itly. Boundary conditions may specify adiabatic walls, specified wall temper- 
ature, or specified heat flux through the walls. 

Cases of practical Interest will have a turbulent layer, of course. For 
this purpose, the Baldwin- Lomax eddy-vlscosl ty model Is used (ref. 4). Pro- 
vision Is made for specifying a constant turbulent Prandtl number for the 
energy equation. 
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Numerical Examples 


Several combinations of geometry, mesh, and flow conditions have been 
calculated for verifying the PARAMAR code. Host of this has been done for 
laminar flow to avoid the possibility of Inaccuracy due to the turbulence 
model. Host of the geometry has been for straight ducts of various cross- 
sections, at varying angles, and for accelerating and decelerating flow. In 
addition some curved ducts and turbomachine blades have been run, which are 
given as examples below. 

Laminar flow In a curved duct . - The duct Is shown In figure 3. Laser 
veloclmeter measurements have been made by T3ylor, Whltelaw, and Ylaneskls 
(ref. 12). The measurements and calculations were made at a Reynolds number 
of 790, based on the hydraulic diameter of the duct. The calculation was 
started at 0.1 m upstream of the start of the bend with uniform flow. The 
mesh used for the calculation by PARAHAR Is shown In figure 4. The Invlscld 
pressure field shown In figure 5 was obtained by using the Denton code 
(ref. 10). Some dlfflcutly was had with negative velocities on the pressure 
surface near the start of the bend. There Is a strong adverse pressure gradi- 
ent In this region, as shown In figure 5, which resulted In negative stream- 
wise velocities being calculated. These negative velocities are closer to the 
wall than the experimental measurements shown In Ref. 12. When the streamwlse 
velocities become negative, the streamwlse momentum Is Increased to a small 
value (2 percent of the freestream value for this case). Thus a small amount 
of momentum Is being artificially added to the boundary layer. The added 
energy Is negligible, but allows calculation through the small reverse flow 
region. This Is a variation on the commonly used Reyhner and Flugge-Lotz 
approximation (ref. 13). 

Carpet plots of streamwlse velocities and cross-section vector plots are 
shown In figures 6 and 7. The results agree qualitatively with the experi- 
mental results. In order to make a detailed comparison with the experiment, 

It Is necessary to Impose the proper flow field at the start of the bend, 
which has not yet been done. 

flow through an axial stator . - The blade passage Is shown In figure 8. 

A cusp has been added to the very blunt leading edge. Flow for this geometry 
nas been calculated for laminar flow at a Reynolds number of 197, based on the 
upstream hydraulic diameter. An Invlscld pressure fie'ld was obtained by the 
ME R l DL and TSONIC codes (refs. 8 and 9). Calculated velocity magnitudes and 
cross-section vector plots are shown In figures 9 and 10. Although experi- 
mental comparisons have not yet been made, the results are qualitatively cor- 
rect. An effort Is being made to obtain a turbulent flow solution at a much 
higher Reynolds number. 
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Figure 3. - Dimensions of curved duct 
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Figure 4. - Computational mesh for curved duct. 
Every fifth marching station shown. Every other 
streamwlse mesh line shown. 



(a) Suction surface. 

Figure 5. - Inviscld pressure for curved duct. 




DISTANCE ALONG WALL, m 
(b) Pressure surface. 
Figure 5. -Concluded. 
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(a) Station 100. (See figure 4. ) 

Figure 6. - Streamwise velocities for curved duct 




Figure 6. - Concluded. 
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PRESSURE WALL 
(a) Station 100. (See figure 4 I 
Figure 7. - Secondary velocity vectors for cruved duct 
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Figure 7. - Continued. 
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PRESSURE WALL 
(d) Station 280. 


Figure 7. - Concluded 



Figure 8. - Axial stator. 
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(c) Station 90. 
Figure?. - Concluded. 
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(a) Station 50. (S? 3 figure & ) 

Figure 10. - Secondary velocity vectors for axial stator. 
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(b) Station 70. 
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(c) Station 90. 
Figure 10. - Concluded. 



